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Abstract 

To study convergence of SMACOF we introduce a modification mSMACOF that 
rotates the configurations from each of the SMACOF iterations to principal components. 

This modification, called mSMACOF, has the same stress values as SMACOF in each 
iteration, but unlike SMACOF it produces a sequence of configurations that properly 
converges to a solution. We show that the modified algorithm can be implemented by 
iterating ordinary SMACOF to convergence, and then rotating the SMACOF solution 
to principal components. The speed of linear convergence of SMACOF and mSMACOF 
is the same, and is equal to the largest eigenvalue of the derivative of the Guttman 
transform, ignoring the trivial unit eigenvalues that result from rotational indeterminacy. 
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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory deleeuwpdx.net/pubfolders/converge has a pdf 
version, the bib hie, the complete Rmd hie with the code chunks, and the R source code. 

1 Introduction 

In (Euclidean, metric, least squares) multidimensional scaling we minimize the real valued 
loss function 

= \ E (1) 

l<i< 7 <n 

defined on M nxp , the space of all n x p matrices. We follow Kruskal (1964) and call cr(X) the 
stress of configuration X. Minimizing stress over p-dimensional configurations is the pMDS 
problem. 

In (1) the matrices of weights W = {w l3 } and dissimilarities A = {//_,} are symmetric, non¬ 
negative, and hollow (zero diagonal). The matrix-valued function D(X) = {dij(X)} contains 
Euclidean distances between the rows of the configuration X, which are the coordinates of n 
points in M p . Thus D(X) is also symmetric, non-negative, and hollow. 

2 Notation 

First some convenient notation, hrst introduced in De Leeuw (1977). Vector e* has n elements, 
with element i equal to +1, and all other elements zero. A V) is the matrix (e^ — ej)(e, — ef)', 
which means elements ( i,i ) and (j,j) are equal to +1, while ( i,j ) and (j, i) are —1. Thus 

d%(X) = (e, - eiYXX'iet - e 3 ) = tr X'A tJ X = tr A tJ C, (2) 


with C = XX'. 

We also define 

V = w ijA-ij, 

and the matrix-valued function B with 


B(X)= E 


& 


Wj 


A • • 


dij(X)> 0 dij(X) 

and B(X) = 0 if X = 0. If we assume, without loss of generality, that 

1 






then 

a(X) = 1 - tr X'B(X)X + Ar X'VX. 


( 3 ) 

( 4 ) 


( 5 ) 
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We also suppose, without loss of generality, that W is irreducible , so that the pMDS problem 
does not separate into a number of smaller pMDS problems. For symmetric matrices 
irreducibity means that we cannot find a permutation matrix II such that IFWII is the direct 
sum of a number of smaller matrices. 

V is symmetric with non-positive off-diagonal elements. It is doubly-centered (rows and 
columns add up to zero) and thus weakly diagonally dominant. It follows that it is positive 
semi-definite (see Varga (1962), section 1.5). Because of irreducibility it has rank n — 1, and 
the vectors in its null space are all proportional to e, the vector with all elements equal to 
+1. The matrix B(X) is also symmetric, positive semi-definite, and doubly-centered for each 
A". It may not be irreducible, because for example B( 0) = 0. 


3 SMACOF 


Define the Guttman Transform of configuration X as 

T(X) = V + B(X)X 


( 6 ) 


The SMACOF algorithm is 

X (k+1) = r(jW) = r fc (X°) (7) 

De Leeuw (1977) shows that the SMACOF iterations (7) tend (in a specifc sense described 
later in this paper) to a fixed point of the Guttman transform, i.e. to an X with T(A) = X. 
If stress is differentiable at X then the fixed points of the Guttman transform are stationary 
points, where the derivative of stress vanishes. Also V + B(X)X = X shows that the columns 
of X are eigenvectors of V + B(X), with eigenvalues equal to one. 

The algorithm was proposed by Guttman (1968), by differentiating stress and setting the 
derivatives equal to zero. This ignores the problem that stress is not differentiable if one or 
more distances are zero, and it also does not say anything about convergence of the algorithm. 
In De Leeuw (1977) a new derivation of the algorithm was given, using a majorization 
argument based on the Cauchy-Schwarz inequality. This make it possible to avoid problems 
with differentiability, and it leads to a simple convergence proof. 


4 Global Convergence of SMACOF 

Following De Leeuw (1977) we also define 

p(X) = £ w i j8ijd ij (X) = tr X'B(X)X, 

l<2<j<n 

V 2 (X) = J2 WijdliX) = tr X'VX; 

l<i<j<n 


and 

MX) 


The main result in De Leeuw (1977) is: 


fit-V) 

ri(X) ■ 


( 8 ) 

( 9 ) 

( 10 ) 
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1. The SMACOF iterates X^ are in the compact set rj(X) < 1, 

2. (p(X fc )}, {//X fc )} and (A(X fc )} are increasing sequences, with limits Poo,Voo and ^oo, 

where = yfpff. 

3. {cr(X fe )} is a decreasing sequence converging to a = 1 — p^. 

4. At any accumulation point X^ of {<x(X fc )} we have a(X 00 ) = and X <*, = T(X 00 ). 

5. p(X ( ' k+1 ' 1 — X^) —* 0 and thus either {(x(X fc )} converges or the set of accumulation 
points of {(r(X fc )} is a continuum. 

So, in words, the sequence of stress values converges monotonically. The sequence of 
configurations is asymptotically regular and has one or more accumulation points. Each 
accumulation point is a fixed point of the Guttman transform. All accumulation points have 
the same function value, which is also the limit of the stress sequence. 

But it is important to emphasize that De Leeuw (1977) did not show that the configuration 
sequence {A^O} actually a Cauchy sequence, and converges to a configuration X^. This is 
basically because of the rotational invariance of the pMDS problem. If K is a rotation matrix, 
i.e. K'K = KK' = /, then dij(XK) = dij(X) for all i and j, and thus cr(XK) = <?(X). If X 
is a fixed point of the Guttman tranform, then so is XK. Consequently there are no isolated 
fixed points, each fixed point is part of a nonlinear continuum of fixed points. 

It is also of interest that De Leeuw (1984) showed that if X is a local minimizer of stress 
then dij(X) > 0 for all i and j such that WijSij > 0. Thus, if weights and dissimilarities are 
positive, stress is differentiable at a local minimum. De Leeuw (1993) shows that stress has 
only a single local maximum at X = 0, and consequently the only possible stationary points 
are local minima and saddle points. It is shown by De Leeuw (2019) that all fixed points of 
the Guttman transform which are not of full column rank are saddle points. Thus if X is a 
local minimizer of pMDS, then adding q zero columns to X produces a saddle point [A" | 0] 
for (p+q)MDS. At saddle points there are always directions in which stress can be decreased, 
and consequently the SMACOF algorithm with enough iterations will eventually get close to 
a continuum of local minima. 


5 Local Convergence of SMACOF 

De Leeuw (1988) studies the local convergence of SMACOF, i.e. the rate of convergence of 
the SMACOF iterations. There is, however, a problem not adequately addressed in that 
article, which has quite a bit of hand-waiving, and some incorrect statements. If there is no 
convergence to a single configuration then the usual rate of convergence is not defined either. 

We know the sequence (? 7 (A^ fc+1 ) — X^)} tends to zero. We can measure its rate of 
convergence by the ratio factor 


and the root factor 


r](X( k+ V - X^) 

Qk ~ rj(XW - ATM) ’ 

(ii) 

r k = r/(X (fc+1) - X {k) ) l/k . 

(12) 
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There is no guarantee that either {q k } or {r k } converges, but 

lim inf q k < lim inf r k < lim sup r k < lim sup q k (13) 

k~ ^oo k— »oo fc—^oo k— >oo 

Thus if linifc^oo q k exists then lim^^ r k exist and the two limit values are equal. See Rudin 
(1976), p. 68, theorem 3.37. 

The rate of convergence of the SMACOF iterations at a fixed point X is computed in De 
Leeuw (1988) as the spectral norm k(X) = ||DT(A") Hoc, i.e. as the modulus of the largest 
eigenvalue of the derivative of the Guttman transform. There are good reasons to compute 
this quantity. From Ostrovski’s Theorem (Ostrowski (1973), theorem 22.1, p. 151 , Ortega 
and Rheinboldt (1970), theorem 10.1.3, p. 300) we know that if n(X) < 1, then X is a point 
of attraction of the SMACOF iteration. If we start close enough then the iterations will 
converge to A". We also know (Ostrowski (1973), theorem 22.2, p. 152) that if n(X) > 1 then 
X is a point of repulsion , which means the iterations will diverge when started in a certain 
solid angle with A" at its apex. 

If we define, with Ortega and Rheinboldt (1970), chapter 9, for a sequence {A"®} converging 
to X, the ratio and root convergence factors 

Qi({v (t) }) = < 14) 

and 

i?i({A (fc) }) = lim sup r){X {k) - X) l/k , (15) 

k—>oo 

then k(X) < 1 implies Ri(X^) = k(X) (Ortega and Rheinboldt (1970), theorem 10.1.4, 
p. 301). Alo Ri(X^) is independent of the norm we have chosen, which is r] in the SMACOF 
case (Ortega and Rheinboldt (1970), theorem 9.9.2, p. 288). Because <5i(A^ fc )) depends on 
the norm we can only assert that Q\(X^) > R\(X^) = k(X), although it is guaranteed in 
the SMACOF case that some norm exists for which there is equality (Ortega and Rheinboldt 
(1970), Notes and Remarks NR 10.1-5, p. 306). 

There is a problem, however, with applying the Ostrowski and Ortega-Rheinboldt results 
in our case. Because of the invariance of distances under rotation we have k(X) = 1, no 
matter what X is. And, basically for the same reason, we have not shown that the SMACOF 
iterations converge to a fixed point at all. De Leeuw (1988) on page 171 acknowledges the 
problem. In fact, he proposes a way around the problem in the proof of theorem 3 on page 175, 
but the proof is not very convincing, although convincing enough to get past the reviewers. I 
will try to be more specific and precise. 

Observe that 

VT{X) = V + V 2 p{X ), (16) 

and 

V 2 a(X) = V — V 2 p(X) = V(I - Vr(X)). (17) 

Because p is convex T> 2 p(X) is symmetric and positive semi-definite. The eigenvalues of 
V r(A) are the generalized eigenvalues of the positive semi-definite matrix pair ( T> 2 p(X ), V). 
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An explicit formula for the derivatives of the Guttman transform, or equivalently for the 
second derivatives of stress, was first given by De Leeuw (1988). We write T>Y(X)\Y] for 
the derivative evaluated at Y. Partial derivatives can be obtained by choosing Y = e,e(,. Of 
course we assume here that djj(X) > 0 for all i < j with WijSij > 0. 

VY{X) [Y] = V + {B(X)Y - if (A, Y)X}, (18) 


with 


H(X,Y) 


Z Wij 


% tr X'AijY 

<ki(X) d%(X) « 


(19) 


These formulas, together with equations (16) and (17), prove the main results in De Leeuw 
(1988). 

1. All eigenvalues of T>Y{X) are non-negative. 

2. (Homogeneity) X is an eigenvector of T>Y(X) with eigenvalue zero. 

3. (Translation) VY(X) has at least p additional eigenvalues equal to zero, corresponding 
with eigenvectors of the form ee' sl where e has all elements equal to one and e s has one 
of its elements equal to one and the others equal to zero. 

4. (Rotation) If X is a fixed point of the Guttman transform then T>T(X) has at least 
\p(j) — 1) eigenvalues equal to one, corresponding with eigenvectors of the form XA, 
with A anti-symmetric. 

5. At a local minimum point X of stress all eigenvalues of DY(X) are less than or equal 
to one. At a strict local minimum they are strictly less than one. 

6. At a saddle point X of stress the largest eigenvalue of T>Y(X) is larger than one. 


6 Modified SMACOF 

To avoid the problems caused by a continuum of fixed points De Leeuw (1988) proposes to 
rotate each SMACOF iterate to a unique position, for example to principal components. For 
any X with singular value decomposition X = KALI we define 11(A") = K A = XL. Note 
that if one or more singular values are equal then II is not unique and becomes a point-to-set 
map. We can rotate within each of the spaces corresponding to multiple singular values, and 
thus we have not completely eliminated indeterminacy due to rotation. In the sequel we 
simply assume that the p singular values of X are different. 

Now consider the iterations 


A (fe+1) = n(r(A (fc) )) = (nr) fc (w°). ( 20 ) 

Let’s call this modified SMACOF, or mSMACOF for short. 

It looks initialy as if mSMACOF require a great deal of extra computation, but this is actually 
not the case. If M is a rotation matrix we have II(AM) = 11(A) and Y(XM) = Y{X)M. 
This implies that Iirn(X) = nr(A) for all X. As a result the sequences {A - ^} and {A^} 
have a simple relationship. If we start the {A^} sequence with X° then for each k we have 
cr(A^) = ct(ATD) and A® = n(A^ fc ^). Thus the two sequences of stress values are exactly 
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the same for SMACOF and mSMACOF and the two configurations, and consequnetly also the 
accumulation points of the two sequences of configurations, just differ by a rotation. To get 
X we can just compute the SMACOF iterate and rotate it to principal components. 

The derivative of Fir is, by the chain rule, 


VUT(X) = VTUT(X))VT(X). (21) 

After some computation we find 

VU(X) \Y] = YL + RAM, (22) 


where X = KAL' and M is the anti-symmetric matrix with off-diagonal elements 


with U = K'YL. 


rriij = 


T XjUji 


(23) 


If A' satisfies n(A) = X then L = 1 and thus M1(X)[Y] = Y + X M and U = K'Y. We 
proceed to compute the eigenvectors and eigenvalues of the Jacobian at an X with II(Af) = X. 

1. If Y = XA with A antisymmetric then M = — L'AL and T>Ii{X)\Y] = 0. This 
corresponds with \p(p — 1) zero eigenvalues 

2. If Y = KA~ l A with A antisymmetric then M = 0. This gives \p(p — 1) eigenvectors 
with eigenvalue one. 

3. If Y = KA~ l D wth D diagonal then M — 0. This gives another p eigenvectors with 
eigenvalue one. 

4. If Y = K ± S with A'j_ a basis for the orthogonal complement of X then M = 0 and we 
have another p(n — p) eigenvalues equal to one. 

This is a complete set of eigenvectors. It turns out all eigenvalues are equal to one, except 
for \p(p — 1) zero eigenvalues. It follows that fixed points of nr are of the form X = XL, 
with X a fixed point of T and L the rotation of X to principal components). Thus T(A") = 
r(X)L = XL = X, and X is a fixed point of both T and 1JT. The eigenvalues of T>TVT(X) 
are the same as the eigenvalues of VY(X), except for the \p{p — 1) eigenvalues equal to one, 
corresponding with XA for antisymmetric A, which are replaced by zero eigenvalues. The 
Ost.rowski and Ortega-Rheinboldt results now apply directly to Fir. If there are no zero 
distances, if the singular values of X are different, and if the largest eigenvalue of JJIir(A") 
is less than one, then X is a point of attraction, the SMACOF iterations converge to X if 
started sufficiently close to it (in the attraction ball of X), and Ri(X) is equal to the largest 
eigenvalue. 


7 Software 

The appendix has an ad-hoc version of the smacof function. Its arguments are 
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args(smacof ) 


## function (delta, w, p = 2, xold = torgerson(delta, p), pea = FALSE, 

## verbose = FALSE, eps = le-15, itmax = 10000) 

## NULL 

The function expects both dissimilarities and weights to be of class dist. The default initial 
configuration uses classical MDS (Torgerson (1958)). We iterate until either the distance 
rj(XW — X^ fc ~b) between successive configurations is less than eps, which has the ridiculously 
small default value of 1CU 15 , or until the maximum number of iterations itmax is reached, 
which defaults to the ridiculously large default value of 10000. If pea is TRUE all iterates 
are rotated to principal components. If verbose is TRUE we print, for each iteration, the 
iteration number k, the stress a(X^), the change rj(X^ — X b), and the root and ratio 
convergence factors. 

The second hie derivative.R computes the Jacobian of the iteration functions at the limit 
(i.e. at the point where SMACOF stops). There are actually six functions in the hie. The 
three functions dGammaA, dPiA and dPiGammaA use the analytical expressions we have derived 
earlier for VT(X), 'DII(X) and Dlir(X). The corresponding functions dGammaN, dPiN and 
dPiGammaN numerically compute the Jacobian using the jacobian function from the numDeriv 
package (Gilbert and Varadhan (2019)). They are basically used to check the formulas. 

8 Examples 

8.1 Ekman 

Time for a numerical example. We will use the color similarity data of Ekman (1954), taken 
from the SMACOF package (De Leeuw and Mair (2009)). We transform Ekman’s similarities 
Sij to dissimilarities 8jj using 8^ = (1 — s^) 3 , and we use unit weights in W. 

SMACOF requires 52 iterations for convergence. The minimum of stress is 0.0110248119, the 
root factor is equal to 0.5099244097 and the ratio factor is 0.6093903884. The eigenvalues of 
VT(X) are 

## 0.999999999999999 
## 0.538510668196408 
## 0.532498554224166 
## 0.529669334191043 
## 0.525541097754551 
## 0.519602718218806 
## 0.516533687841055 
## 0.513727908691609 
## 0.510295310409166 
## 0.506357695934493 
## 0.495649713446001 
## 0.481518780073304 


## 0.469548625536452 
## 0.445038589020916 
## 0.410479252654483 
## 0.394284315478173 
## 0.357670750114585 
## 0.348395413045201 
## 0.330341477528788 
## 0.313520384427864 
## 0.287598015841956 
## 0.280399104121623 
## 0.267278474596759 
## 0.247631495462991 
## 0.216576009083046 
## 0.000000000000000 
## 0.000000000000000 
## 0.000000000000000 

As an aside, the Ekman example (with this particular transformation of the similarities) 
is special, because the two-dimensional solution is actually the global minimum over all 
configurations (i.e. the global mnimum of pMDS for all 1 < p < n). As shown in De Leeuw 
(2014) this follows from the fact that the two unit eigenvalues of V + B(X) are actually its two 
largest eigenvalues, and consequently A" is the solution of the convex full-dimensional nMDS 
relaxation of the pMDS problem (De Leeuw, Groenen, and Mair (2016)). The eigenvalues of 
V+B(X) are 

## 1.000000000000000 
## 0.999999999999999 
## 0.923497086367287 
## 0.907901212970889 
## 0.862936584878895 
## 0.852692003103344 
## 0.829803620857356 
## 0.814556167662381 
## 0.793238576338501 
## 0.791651722456648 
## 0.786442678118770 
## 0.747679475705025 
## 0.728268247434340 
## 0.000000000000001 

We now set pca=TRUE and run mSMACOF. It requires 52 iterations for convergence. The 
minimum of stress is 0.0110248119, and the root factor is equal to 0.5099406626 and the 
ratio factor is 0.6146787643. The eigenvalues of VHY^X) at the fixed point X are 

## 0.538510668196406 
## 0.532498554224166 
## 0.529669334191042 
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## 0.525541097754553 
## 0.519602718218806 
## 0.516533687841055 
## 0.513727908691609 
## 0.510295310409167 
## 0.506357695934493 
## 0.495649713446000 
## 0.481518780073304 
## 0.469548625536451 
## 0.445038589020917 
## 0.410479252654483 
## 0.394284315478172 
## 0.357670750114584 
## 0.348395413045201 
## 0.330341477528788 
## 0.313520384427864 
## 0.287598015841956 
## 0.280399104121622 
## 0.267278474596759 
## 0.247631495462990 
## 0.216576009083047 
## 0.000000000000001 
## 0.000000000000001 
## 0.000000000000000 
## 0.000000000000000 

They are equal to the eigenvalues of T>T{X) and the root convergence factor ^({AT^}) is 
0.5385106682. 

8.2 De Gruijter 

The second example are dissimilarties between nine Dutch political parties, collected a long 
time ago by De Gruijter (1967). 


## 


KVP 

PvdA 

VVD 

ARP 

CHU 

CPN 

PSP BP 

## 

PvdA 

2.63 







## 

VVD 

2.27 

3.72 






## 

ARP 

1.60 

2.64 

2.46 





## 

CHU 

1.80 

3.22 

1.97 

0.20 




## 

CPN 

4.54 

2.12 

5.13 

4.84 

4.80 



## 

PSP 

3.73 

1.59 

4.55 

3.73 

4.08 

1.08 


## 

BP 

4.18 

4.22 

3.90 

4.28 

3.96 

3.34 

3.88 

## 

D66 

3.17 

2.47 

1.67 

3.13 

3.04 

4.42 

3.36 4.36 


We analyze the data with SMACOF using three dimensions. 

SMACOF requires 779 iterations for convergence. The minimum of stress is 0.003442194, the 
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root factor is equal to 0.9565805061, and the ratio factor is 0.9580406219. The eigenvalues of 
VT(X) are 

## 1.000000000000000 
## 1.000000000000000 
## 1.000000000000000 
## 0.965505429805660 
## 0.940592046981168 
## 0.919047686446446 
## 0.863993920924432 
## 0.822263696226350 
## 0.810020971414307 
## 0.771947328045072 
## 0.728539213663602 
## 0.712621684571534 
## 0.659318624263221 
## 0.638485365688917 
## 0.624552464148481 
## 0.616085432872672 
## 0.557480497708779 
## 0.501954186928525 
## 0.458630667850016 
## 0.450598367846592 
## 0.355243400669833 
## 0.309186479174060 
## 0.247708397109092 
## 0.000000000000002 
## 0.000000000000001 
## 0.000000000000000 
## 0.000000000000000 

In this case there is no guarantee that we have the three-dimensional global minimum. The 
unit eigenvalues of the matrix V + B(X) are not the three largest ones. 

## 1.079524009371954 
## 1.032606649163671 
## 1.000000000000000 
## 1.000000000000000 
## 1.000000000000000 
## 0.986706272372898 
## 0.971839080877692 
## 0.906211919383163 
## 0.000000000000002 

mSMACOF in three dimensions requires 783 iterations for convergence. The minimum 
of stress is 0.003442194, the root factor is equal to 0.9568196354, and the ratio factor is 
0.8713099253. The eigenvalues of VHT^X) are 


11 



## 0.965505429805661 
## 0.940592046981166 
## 0.919047686446444 
## 0.863993920924432 
## 0.822263696226350 
## 0.810020971414307 
## 0.771947328045072 
## 0.728539213663603 
## 0.712621684571535 
## 0.659318624263221 
## 0.638485365688918 
## 0.624552464148482 
## 0.616085432872673 
## 0.557480497708779 
## 0.501954186928524 
## 0.458630667850017 
## 0.450598367846592 
## 0.355243400669832 
## 0.309186479174062 
## 0.247708397109091 
## 0.000000000000001 
## 0.000000000000001 
## 0.000000000000001 
## 0.000000000000000 
## 0.000000000000000 
## 0.000000000000000 
## 0.000000000000000 

Again, as expected, the eigenvalues for SMACOF and mSMACOF are the same and the root 
convergence factor is the spectral norm of the mSMACOF derivative, which is 0.9655054298. 


9 Conclusion 

From a practical point of view there is no need to ever use modified SMACOF. We only use 
the Guttman transform, and compute the largest eigenvalue of its derivative at the solution 
X, ignoring the \p(p — 1) trivial unit eigenvalues. That largest eigenvalue, if it is strictly 
smaller than one, is the root convergence factor. In that case SMACOF can be said to 
converge to II(A"), which is the SMACOF solution rotated to principal components. 


12 



10 Appendix: Code 

10.1 smacof.R 


torgerson <- function (delta, p = 2) { 
h <- as .matrix(delta ~ 2) 
n <- nrow (h) 
j <- diag (n) - (1 / n) 

h <- -(j y.*y. h y.*y. j) / 2 

e <- eigen (h) 

return (e$vectors[, 1:p] diag(sqrt (e$values [1 :p]))) 

} 

smacof <- 

function (delta, 
w, 

p = 2, 

xold = torgerson (delta, p), 

pea = FALSE, 

verbose = FALSE, 

eps = le-15, 

itmax = 10000) { 

n <- round ((1 + sqrt (1 + 8 * length (delta))) / 2) 

v <- -as.matrix (w) 

diag(v) <- -rowSums(v) 

vinv <- solve(v + (1 / n)) - (1 / n) 

delta <- delta / sqrt (sum (w * delta 2) / 2) 

dold <- dist (xold) 

sold <- sum (w * (delta - dold) 2) / 2 
eold <- Inf 
itel <- 1 
repeat { 

b <- as.matrix (-w * delta / dold) 
diag (b) <- -rowSums(b) 
xnew <- vinv °/o*°/ 0 b %*% xold 
if (pea) { 

xsvd <- svd (xnew) 
xnew <- xnew xsvd$v 

} 

dnew <- dist (xnew) 

snew <- sum (w * (delta - dnew) 2) / 2 
enew <- sqrt (sum (v * tcrossprod(xold - xnew))) 
rnew <- enew (1 / itel) 
qnew <- enew / eold 
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if (verbose) { 
cat ( 

"itel ", 

f ormatC(itel , digits = 4, format = 
"loss ", 

f ormatC(snew, digits = 15, format = 
"chan ", 

formate (enew, digits = 15, format = 
"renf ", 

f ormatC(rnew, digits = 15, format = 
"qcnf ", 

f ormatC(qnew, digits = 15, format = 
"\n" 


"d"), 
"f ") 
"f ") 
"f ") 
"f ") 


) 

> 

if ((enew < eps) I I (itel == itmax)) 
break 

xold <- xnew 
dold <- dnew 
sold <- snew 
eold <- enew 
itel <- itel + 1 


> 


out <- 
list ( 

itel = itel, 
x = xnew, 
s = snew, 
q = qnew, 
r = rnew, 
b = vinv %*% b 

) 

return (out) 

} 


10.2 derivative.R 

dGammaA <- function (x, delta, w) { 
n <- nrow (x) 
p <- ncol (x) 

h <- matrix (0, n * p, n * p) 
d <- dist (x) 

delta <- delta / sqrt (sum (w * delta 2) / 2) 
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dmat <- as.matrix (d) 
wmat <- as.matrix (w) 
emat <- as.matrix (delta) 
v <- -wmat 

diag(v) <- -rowSums (v) 
b <- as.matrix(-w * delta / d) 
diag(b) <- -rowSums (b) 
vinv <- solve (v +(l/n))-(l/n) 
for (s in l:p) { 
for (t in l:p) { 

gst <- matrix (0, n, n) 
for (i in l:n) { 
for (i in l:n) { 
if (i == j) 
next 

gs <— x [i, s] - x[j, s] 
gt <- x[i, t] - x[j, t] 
gst[i, j] <- 

-wmat[i, j] * emat[i, j] * gs * gt / (dmat[i, j] 3) 

} 

} 

diag(gst) <- -rowSums (gst) 

h[(s - 1) * n + l:n, (t - 1) *n + l:n] <- -vinv 0 / 0 * 0 / 0 gst 

> 

} 

for (s in 1: p) { 

kn<- (s - 1) *n+ 1 :n 

h[kn, kn] <- h[kn, kn] + vinv °/„*°/„ b 

> 

return (h) 

} 

dPiA <- function (x) { 
n <- nrow (x) 
p <- ncol (x) 
sx <- svd (x) 
xu <- sx$u 
xv <- sx$v 
xd <- sx$d 

h <- matrix (0, n * p, n * p) 
for (s in l:p) { 
for (i in l:n) { 

ir <- (s - 1) * n + i 
e <- matrix (0, n, p) 
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m <- matrix (0, p, p) 
e[i, s] <- 1 

u <- crossprod (xu, e %*% xv) 
for (k in 1: p) { 
for (1 in 1: p) { 
if (k == 1) 
next 

m[k, 1] <- 

(xd[k] * u[k, 1] + xd[l] * u[l, k]) / (xd[k] “ 2 - xd[l] “ 2) 

} 

} 

h[, ir] <- as.vector (e %*% xv - xu °/ 0 *°/ 0 diag (xd) %*% m) 

> 

} 

return (h) 


dGammaN <- function (x, delta, w) { 
n <- nrow (x) 
p <- ncol (x) 

delta <- delta / sqrt (sum (w * delta 2) / 2) 
v <- - as.matrix (w) 
diag (v) <- - rowSums(v) 
vinv <- solve (v + (1 / n)) - (1 / n) 
guttman <- function (x) { 
d <- dist (matrix (x, n, p)) 
b <- - as.matrix (w * delta / d) 
diag (b) <- - rowSums (b) 
z <- vinv %*% b %*% matrix (x, n, p) 
return (as.vector (z)) 

} 

return (jacobian (guttman, as.vector (x))) 

} 

dPiGammaN <- function (x, delta, w) { 
n <- nrow (x) 
p <- ncol (x) 

delta <- delta / sqrt (sum (w * delta 2) / 2) 
v <- - as.matrix (w) 
diag (v) <- - rowSums (v) 
vinv <- solve (v + (1 / n)) - (1 / n) 
guttman <- function (x) { 
d <- dist (matrix (x, n, p)) 
b <- - as.matrix (w * delta / d) 
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diag (b) <- - rowSums (b) 
z <- vinv 1*1 b 1*1 matrix (x, n, p) 
z <- z 1*1 svd(z)$v 
return (as.vector (z)) 

} 

return (jacobian (guttman, as.vector (x))) 

} 

dPiGammaA <- function (x, delta, w) { 
n <- nrow (x) 

guttman <- function (x, delta, w) { 
d <- dist (x) 

b <- as.matrix (- w * delta / d) 

diag (b) <- - rowSums (b) 

v <- - as.matrix (w) 

diag (v) <- - rowSums (v) 

vinv <- solve(v + (1 / n)) - (1 / n) 

return (vinv 0 / 0 *°/ 0 b 1*1 x) 

> 

h <- dGammaA (x, delta, w) 
g <- dPiA (guttman (x, delta, w)) 
return (g 1*1 h) 

} 
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